%pip install ydata_profiling numpy pandas matplotlib seaborn scikit-learn PyQt6 ipykernel

%matplotlib inline

import os
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from ydata_profiling import ProfileReport
from sklearn.preprocessing import LabelEncoder

data = pd.read_csv('weatherAUS.csv')

report_file = 'weatherAUS_report.html'

if not os.path.exists(report_file):
    profile = ProfileReport(data, title = 'Weather in Australia - Report')

    profile.to_file('weatherAUS_report_1.html')
    profile.dump('report_1')

profile = ProfileReport().load('report_1.pp')
profile.to_notebook_iframe()

numeric_cols = data.select_dtypes(include = ['float64']).columns

plt.figure(figsize=(25, 15))  

for column in range(len(numeric_cols)):
    plt.subplot(4, 4, column + 1)
    sns.boxplot(y = data[numeric_cols[column]])  
    plt.title(f'The boxplot of {numeric_cols[column]}')  
    
plt.tight_layout()  
plt.show()

for column in numeric_cols:

    Q1 = data[column].quantile(0.25)
    Q3 = data[column].quantile(0.75)
    mean = data[column].median()

    IQR = Q3 - Q1

    upper = Q3 + IQR * 1.5
    lower = Q1 - IQR * 1.5

    data.loc[data[column] >= upper, column] = np.nan
    data.loc[data[column] <= lower, column] = np.nan

plt.figure(figsize=(25, 15))  

for column in range(len(numeric_cols)):
    plt.subplot(4, 4, column + 1)
    sns.boxplot(y = data[numeric_cols[column]])  
    plt.title(f'The boxplot of {numeric_cols[column]}')  
    
plt.tight_layout()
plt.show()

data.sort_values(by = 'Date', inplace = True, ascending = True, ignore_index = True)

dates = pd.to_datetime(data['Date'])

data.drop('Date', axis = 1, inplace = True)
data.dropna(subset = 'RainTomorrow', inplace = True)

data.insert(0, 'Day', dates.dt.day.astype('float64'))
data.insert(0, 'Year', dates.dt.year.astype('float64'))
data.insert(0, 'Month', dates.dt.month.astype('float64'))

data['RainToday'] = data['RainToday'].map({'Yes': 1, 'No': 0}).astype('float64')
data['RainTomorrow'] = data['RainTomorrow'].map({'Yes': 1, 'No': 0}).astype('float64')

data.head()

leng = data.__len__()

data_train_set = data.iloc[: np.ceil(leng * 0.65).astype(int)].copy()
data_val_set = data.iloc[np.ceil(leng * 0.65).astype(int) : np.ceil(leng * 0.8).astype(int)].copy()
data_test_set = data.iloc[np.ceil(leng * 0.8).astype(int) :].copy()

le = LabelEncoder()

for feature in data_train_set.select_dtypes('object').columns.to_list():
    data_train_set.loc[:, feature] = le.fit_transform(data_train_set[feature].astype(str))
    data_test_set.loc[:, feature] = le.transform(data_test_set[feature].astype(str))
    data_val_set.loc[:, feature] = le.transform(data_val_set[feature].astype(str))

data_train_set.head()

correlation_matrix = data_train_set.corr()

plt.figure(figsize = (20, 10))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f')
plt.title('Original Correlation Matrix')
plt.show()

correlation_with_target = correlation_matrix['RainTomorrow'].sort_values(ascending = True)

# Identify features with low correlation to the target (absolute correlation < 0.1)
low_correlation_features = correlation_with_target[correlation_with_target.abs() < 0.1].index.tolist()

plt.figure(figsize = (20, 10))
sns.heatmap(correlation_matrix.loc[low_correlation_features], annot=True, cmap='coolwarm', fmt='.2f')
plt.title('Correlation Matrix of Low Correlated Features')
plt.show()

# Drop low correlated features from the dataset
data_train_set.drop(columns = low_correlation_features, inplace = True)
data_test_set.drop(columns = low_correlation_features, inplace = True)
data_val_set.drop(columns = low_correlation_features, inplace = True)

correlation_matrix_dropped = data_train_set.corr()

plt.figure(figsize = (20, 10))
sns.heatmap(correlation_matrix_dropped, annot=True, cmap='coolwarm', fmt='.2f')
plt.title('Updated Correlation Matrix after Dropping Low Correlated Features')
plt.show()

plt.figure(figsize=(14,5))
missing = data_train_set.isnull().sum() / data_train_set.shape[0] *100
missing = missing.sort_values().plot.bar()

plt.xticks(rotation = 45)
plt.title('Bar plot of missing values with percentage')
plt.show()

%pip install lightgbm

import lightgbm as lgb

def determine_objective(feature):
    if feature.dtype.kind == 'f':
        return 'regression'
    elif len(feature.value_counts()) == 2:
        return 'binary'
    else:
        return 'multiclass'

inputation_models = {}

categorical_features =  data_train_set.select_dtypes(include = ['object']).columns.tolist()

def Training_Inputation_Models(dataset):
    columns = dataset.columns.tolist()
    columns.remove('RainTomorrow')

    for target in columns:
        
        indices = dataset.index[dataset[target].isna()]

        if indices.empty:
            continue

        """
        Dropping columns which target contains NaN values and saving corresponding features in X_target.
        It will be used for predicting dropped NaN values.
        """
        if not indices.empty:
            X_target = dataset.loc[indices].drop(target, axis = 1)

        dropped_data = dataset.dropna(subset = target)

        X_train, Y_train = dropped_data.drop(target, axis = 1), dropped_data[target]

        """
        Parameters of LightGBM model
        """
        objective = determine_objective(dataset[target])

        params = {
            'objective': objective,                  # using regression, binary or multiclass model
            'boosting': 'gbdt',                      # model uses Gradient Boosting which aimed to minimize loss (error)
            'learning_rate': 0.1,                    # the step which model makes in each iteration for minimizing loss
            'num_leaves': 100,                       # number of leaves for tree model
            'verbose': -1,                           # force model not to show additional text information
            'device_type': 'cpu'                     # it says to algorithm to use gpu
        }

        if objective == 'multiclass':
            params['num_class'] = len(dataset[target].unique())

        """
        As there is no test case, we train our model with full dataset (after dropping NaN target values).
        """
        train_set = lgb.Dataset(X_train, label = Y_train)

        model = lgb.train(params, train_set, num_boost_round = 100)
        if not indices.empty:
            Y_target = model.predict(X_target, num_iteration = model.best_iteration)
            
            """
            Writing predicted targets back to the data which also helps to have richer dataset in each iteration.
            """
            dataset.loc[indices, target] = Y_target
        
        inputation_models[target] = model

def Inputation_Of_Missing_Values(dataset):

    nan_contained_cols = dataset.columns[dataset.isna().any()].tolist()

    for target in nan_contained_cols:
        indices = dataset.index[dataset[target].isna()]

        X_target = dataset.loc[indices].drop(target, axis = 1)
        Y_target = inputation_models[target].predict(X_target)

        dataset.loc[indices, target] = Y_target


Training_Inputation_Models(data_train_set)
Inputation_Of_Missing_Values(data_val_set)
Inputation_Of_Missing_Values(data_test_set)

%pip install xgboost

import numpy as np
import matplotlib.pyplot as plt
from xgboost import XGBClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.neural_network import MLPClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, VotingClassifier
from sklearn.model_selection import GridSearchCV, TimeSeriesSplit, RandomizedSearchCV
from sklearn.metrics import f1_score, recall_score, precision_score, precision_recall_curve

X_train, Y_train  = data_train_set.drop('RainTomorrow', axis = 1), data_train_set['RainTomorrow']
X_test, Y_test  = data_test_set.drop('RainTomorrow', axis = 1), data_test_set['RainTomorrow']
X_val, Y_val  = data_val_set.drop('RainTomorrow', axis = 1), data_val_set['RainTomorrow']

tscv = TimeSeriesSplit(n_splits = 5)

def Threshold_Optimizer(model, x_test, y_test):
    y_val_prob = model.predict_proba(X_val)[:, 1]

    precisions, recalls, thresholds = precision_recall_curve(Y_val, y_val_prob)
    f1_scores = 2 * (precisions * recalls) / (precisions + recalls)
    optimal_index = np.argmax(f1_scores)
    optimal_threshold = thresholds[optimal_index]

    """
    Precision-Recall curve illustrates the optimal threshold value which maximizes both recall and precision (as a result f1 score). As long as we do not have special requirement such as bounded recall or precision. We keep this threshold value.
    """

    plt.figure()
    plt.plot(recalls, precisions, color='b',  label='PR curve')
    plt.scatter(recalls[optimal_index], precisions[optimal_index], marker='o', color='black', label='Optimal threshold')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.title('Precision-Recall curve')
    plt.legend(loc="lower left")
    plt.show()

    y_test_prob = model.predict_proba(x_test)[:, 1]

    y_test_pred = (y_test_prob >= optimal_threshold).astype(int)

    print('F1 Test: ', f1_score(y_test, y_test_pred))
    print('Recall Test: ', recall_score(y_test, y_test_pred))
    print('Precision Test: ', precision_score(y_test, y_test_pred))

def Results_Before_Threshold_Optimization(model, x_test, y_test):
    pred_train_set = model.predict(X_train)
    pred_test_set = model.predict(x_test)

    print('Recall Train: ', recall_score(Y_train, pred_train_set))
    print('Recall Test: ', recall_score(y_test, pred_test_set))
    print('Precision Train: ', precision_score(Y_train, pred_train_set))
    print('Precision Test: ', precision_score(y_test, pred_test_set))
    print('F1 Train: ', f1_score(Y_train, pred_train_set))
    print('F1 Test: ', f1_score(y_test, pred_test_set))

pipeline = make_pipeline(
    StandardScaler(),
    KNeighborsClassifier(n_neighbors = 11)
)

grid_search_cv = GridSearchCV(
    pipeline,
    param_grid = {},
    cv = tscv,
    scoring = 'f1',
    n_jobs = -1
)

grid_search_cv.fit(X_train, Y_train)

knn_model = grid_search_cv.best_estimator_

Results_Before_Threshold_Optimization(knn_model, X_val, Y_val)

Threshold_Optimizer(knn_model, X_val, Y_val)

"""
Here we use all possible combinations, while some pentalties do not support all solvers, we should define permutations of parameters such that

params_log = [
    {'log_reg__penalty': ['l1'], 'log_reg__C': [10 ** i for i in range(-3, 4)], 'log_reg__class_weight': ['balanced', None], 'log_reg__solver': ['liblinear', 'saga']},
    {'log_reg__penalty': ['l2'], 'log_reg__C': [10 ** i for i in range(-3, 4)], 'log_reg__class_weight': ['balanced', None], 'log_reg__solver': ['lbfgs', 'liblinear', 'newton-cg', 'newton-cholesky', 'sag', 'saga']},
    {'log_reg__penalty': ['elasticnet'], 'log_reg__C': [10 ** i for i in range(-3, 4)], 'log_reg__class_weight': ['balanced', None], 'log_reg__solver': ['saga']},
    {'log_reg__penalty': [None], 'log_reg__C': [10 ** i for i in range(-3, 4)], 'log_reg__class_weight': ['balanced', None], 'log_reg__solver': ['lbfgs', 'newton-cg', 'newton-cholesky', 'sag', 'saga']}
]

"""

optimal_params = {
    'C': 2,
    'class_weight': 'balanced',
    'penalty': 'l2',
    'solver': 'newton-cholesky'
}

pipeline = make_pipeline(
    StandardScaler(),
    LogisticRegression(**optimal_params)
)

grid_search_cv = GridSearchCV(
    pipeline,
    param_grid = {},
    cv = tscv,
    n_jobs = -1,
    scoring='f1'
)

grid_search_cv.fit(X_train, Y_train)

log_model = grid_search_cv.best_estimator_

Results_Before_Threshold_Optimization(log_model, X_val, Y_val)

Threshold_Optimizer(log_model, X_val, Y_val)

random_params = {
    'n_estimators': 150,
    'max_features': 'sqrt',
    'max_depth': 13,
    'min_samples_split': 2,
    'min_samples_leaf': 4,
    'bootstrap': True,
    'class_weight': 'balanced'
}

grid_search_cv = GridSearchCV(
    RandomForestClassifier(**random_params),
    param_grid = {},
    cv = tscv,
    n_jobs = -1,
    scoring = 'f1'
)

grid_search_cv.fit(X_train, Y_train)

rf_model = grid_search_cv.best_estimator_

Results_Before_Threshold_Optimization(rf_model, X_val, Y_val)

Threshold_Optimizer(rf_model, X_val, Y_val)

xgb_params = {
    'max_depth': 8,
    'n_estimators': 200,
    'learning_rate': 0.072,
    'scale_pos_weight': 2.5,
}

grid_search_cv = GridSearchCV(
    XGBClassifier(**xgb_params),
    param_grid = {},
    cv = tscv,
    n_jobs = -1,
    scoring = 'f1'
)

grid_search_cv.fit(X_train, Y_train)

xgb_model = grid_search_cv.best_estimator_

Results_Before_Threshold_Optimization(xgb_model, X_val, Y_val)

Threshold_Optimizer(xgb_model, X_val, Y_val)

nn_params = {
    'hidden_layer_sizes': (50, 30),
    'activation':  'tanh',
    'solver': 'sgd',
    'learning_rate': 'adaptive',
    'learning_rate_init': 0.001,
    'alpha': 0.01,
    'momentum': 0.99,
    'early_stopping': True,
    'max_iter': 500
}

pipeline = make_pipeline(
    StandardScaler(),
    MLPClassifier(**nn_params)
)

random_search_cv = GridSearchCV(
    pipeline,
    param_grid = {},
    scoring = 'f1',
    cv = tscv,
    n_jobs = -1
)

random_search_cv.fit(X_train, Y_train)

nn_model = random_search_cv.best_estimator_

Results_Before_Threshold_Optimization(nn_model, X_val, Y_val)

Threshold_Optimizer(nn_model, X_val, Y_val)

def Voting_Classifier_Model(x_test, y_test):
    voting_classifier = VotingClassifier(estimators = [
        ('rf', rf_model),
        ('nn', nn_model),
        #('xgb', xgb_model), # Uncomment if you want to include XGBoost model
    ],  voting = 'soft'
    )

    voting_classifier.fit(X_train, Y_train)

    Results_Before_Threshold_Optimization(voting_classifier, x_test, y_test)

    Threshold_Optimizer(voting_classifier, x_test, y_test)

Voting_Classifier_Model(X_val, Y_val)

Voting_Classifier_Model(X_test, Y_test)